function y = MeasureTemp(path)

if exist('path','var')
    cd(path)
else
    path = cd;
end

numRuns = length(dir('*.h5'));

if numRuns == 1
    error('Only one measurement in folder, can''t measure temp.')
elseif numRuns == 0
    error('No measurements in folder')
end

cons = Constants();

for n = 1:numRuns-1
    %     n
    img_time = h5read([path '\run_' num2str(n) '.h5'],'/Image/exposuretime');
    img_var = h5read([path '\run_' num2str(n) '.h5'],'/Image/configuration');
    expp = ParamImgSys(GetImgConfig(img_var));
    % Get variable from Cicero
    t(n) = 1e-3*double(h5read([path '\run_' num2str(n) '.h5'],'/VariableValues/t'));
    
    yC = 1:1040;
    xC = 1:1392;
    
    [pX,pY] = ProcessImage([path '\run_' num2str(n) '.h5'],1,[1 1392 1 1040],1);
    
    % Estimate number of atoms using camera
    amp = (pX(1) + pY(1))/2; % Take mean amplitude of two Gaussian fits.
    %     capture_fraction = 0.5*(1 - sqrt(1 - (expp.na_telec)^2)); % fraction of total light captured by lens ( ~ 0.5 %)
    total_num_electrons = amp*sqrt(2*pi)*mean([pX(3) pY(3)]);
    total_photons = total_num_electrons/(expp.coll_eff);
    num_atoms(n) = 2*total_photons/(img_time*cons.Gamma); % num_atoms = (num_photons/time) / (scattering rate)
    
    % Give the width as the std deviation in um
    px_size = expp.px_size; % um
    mag = expp.mag; % magnification
    wY(n) = pY(3)*px_size/mag; % Std dev in um
    wX(n) = pX(3)*px_size/mag;
    
    posX(n) = pX(2)*px_size/mag;
    posY(n) = pY(2)*px_size/mag;
    
end

tt = t.^2;

wYm2 = 2*(wY*1e-0).^2; % m^2, converted from um.
wXm2 = 2*(wX*1e-0).^2;

coeffsY = polyfit(tt(1:end),wYm2(1:end),1);
coeffsX = polyfit(tt(1:end),wXm2(1:end),1);

v02y = coeffsY(1);
v02x = coeffsX(1);

m = 86.9*1.66e-27; % mass of Rb87
k = 1.38e-23;

Tx = v02x*m/(2*k) % K
Ty = v02y*m/(2*k) % K

r0y = sqrt(coeffsY(2));
r0x = sqrt(coeffsX(2));

figure(2)
plot(tt,wYm2,'.',tt,wXm2,'x',tt,polyval(coeffsX,tt),tt,polyval(coeffsY,tt))
legend('Horizontal','Vertical')

% x = 1/2 * g * t^2
gCX = polyfit(t,posX,2);
g = 2*gCX(1)

tHiRes = linspace(min(t),max(t));

figure(3)
plot(t,posY,'.',t,posX,'x',tHiRes,polyval(gCX,tHiRes))